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ABSTRACT 

We investigate the effects of weakly-interacting massive particle (WIMP) dark matter annihilation on 

the formation of Population III.l stars, which are theorized to form from the collapse of gas cores at 

the centers of dark matter minihalos. We consider the relative importance of cooling due to baryonic 

radiative processes and heating due to WIMP annihilation. We analyze the dark matter and gas profiles 

of several halos formed in cosmological-scale numerical simulations. The heating rate depends sensitively 

^S^ , on the dark matter density profile, which we approximate with a power law p^ oc r~°'^ , in the numerically 

Mh' unresolved inner regions of the halo. If we assume a self-similar structure so that a^ ~ 1.5 as measured 

on the resolved scales ~ 1 pc, then for a fiducial WIMP mass of 100 GeV, the heating rate is typically 

^ ■ much smaller (< 10^'^) than the cooling rate for densities up to nn = 10^^ cm^'^. In one case, where 

"o^ ' a-^ — 1.65, the heating rate becomes similar to the cooling rate by a density of nn = 10^^ cm~^. The 

d . dark matter density profile is expected to steepen in the central baryon-dominated region ^ 1 pc due to 

adiabatic contraction, and we observe this effect (though with relatively low resolution) in our numerical 

models. From these we estimate a^ ~ 2.0. The heating now dominates cooling above njj — 10^** cm""^, 

►^ . in agreement with the previous study of Spolyar, Freese & Gondolo. We expect this leads to the 

Qs ' formation of an equilibrium structure with a baryonic and dark matter density distribution exhibiting a 

\^ I flattened central core. Examining such equilibria, we find total luminosities due to WIMP annihilation 

t^^ i are relatively constant and ^ 10^ Lq, set by the radiative luminosity of the baryonic core. We discuss 

CO ' the implications for Pop III.l star formation, particularly the subsequent growth and evolution of the 

protostar. Even if the initial protostar fails to accumulate any additional dark matter, its contraction to 

}~^ • the main sequence could be significantly delayed by WIMP annihilation heating, potentially raising the 

QQ I mass scale of Pop III.l stars to masses ^ 100 AIq. 

f~^ , Subject headings: cosmology: theory — dark matter — early universe — galaxies: formation — stars: 

formation 



1. INTRODUCTION 



Population III stars are defined to be those whose formation and evolution are independent of metallicity (McKee 
J^ ■ & Tan 2008; O'Shea et al. 2008), since their metallicity is extremely low: close to or equal to that arising from big 
bang nucleosynthesis. Population III.l stars are defined as having their formation be independent of other stars or other 
astrophysical objects, so that their initial conditions are determined solely by cosmology. These stars will be the first 
objects to form in a given region of the universe and they are likely to play an important role in cosmic reionization and 
in laying the foundations for the formation and build-up of galaxies. It is possible that they are the direct or indirect 
progenitors of supermassive black holes. 

Within the commonly accepted ACDM framework. Pop III.l stars form within dark matter halos. Indeed, in those 
halos that do form stars, only one star appears to form in the center of each halo (Abel et al. 2002; Bromm, Coppi, & 
Larson 2002). 

One of the most theoretically well-motivated cold dark matter candidates is a Weakly-Interacting Massive Particle 
(WIMP). Supersymmctric theories with R-parity naturally provide a stable dark matter candidate which could compose 
all or part of the dark matter in the Universe. 

It has been pointed out that, if dark matter consists of a WIMP such as the supersymmetric neutralino, the energy 
released by the annihilation of these particles could influence early structure formation, star formation and protostellar 
evolution (Ripamonti, Mapelli & Ferrara 2007; Ascasibar 2007; Spolyar, Freese, & Gondolo 2008; locco 2008; Freese, 
Spolyar & Aguirre 2008; Freese et al. 2008b, c). The effects on stellar evolution at fixed mass have also been investigated 
(Taoso et al. 2008; Yoon, locco & Akiyama 2008). Spolyar et al. (2008) show that, for the adiabatically-contracted 
Navarro, Frenk & White (1996) (NEW) dark matter density profiles they considered, dark matter heating can overwhelm 
gas cooling in the innermost region of a star-forming minihalo, and they propose that this can then lead to a dark matter 
powered star. 
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In this paper we revisit the scenario investigated by Spolyar et al. (2008). In 121 we derive an analytic expression for 
the dark matter heating rate, including a simplified treatment of radiative transport. In Sj3]we present our results of the 
assessment of the importance of WIMP annihilation heating for several halos formed in numerical simulations of cosmic 
structure formation (O'Shea & Norman 2007). We describe the dark matter density structure in H3.1l and the properties 
of the baryons in H3.'2[ We compare the WIMP annihilation heating rates and the baryonic cooling rates and discuss the 
equilibrium structure of dark matter powered protostars in ^'S.'6\ We discuss the implications for subsequent protostellar 
evolution in ij3.4l We conclude in 311 



2. ANALYTIC PRESCRIPTION FOR WIMP ANNIHILATION HEATING 



Let US consider WIMPs with a mass 



Let < aaV > he the WIMP annihilation cross section times the relative 



velocity, averaged over the momentum distribution of the WIMPs. For numerical evaluations in this paper we will assume 
it has a value 3 x lO^^^cm^ s~^ (e.g. Jungman, Kamionkowski & Griest 1996). The number of photons produced by WIMP 
annihilation per unit volume, per unit time and per unit energy, at the point r' is (Hall & Gondolo 2006; Bergstrom, UUio 
& Buckley 2006): 
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where Px{f') is the WIMP density at r' and dN^/dE^ is the number of photons produced per unit energy per annihilation. 
Let us consider electron scattering at a location r, by photons produced at the location r' . Some of the photons produced 
at r' will scatter off electrons before reaching the location r. The number of photons reaching r without scattering is 
smaller than the number produced at r' by the factor S: 

'a^-^{E^) Jq^ ds ne{s) when r' < r 

cxp — (Tg — f{Ej) J , ^ ds ne(s) ) when r' > r 

where s± = r' cos 6 ± v r^ — r'^ sin^ 9 and (t^-^ is the cross section. Eq. [2] takes account of the fact that photons can pass 
through the sphere r =constant twice when r' > r. We note that S* is a function of the angle 6 between the path of the 
photon and the line joining r' with the center of the spherical cloud. 

Consider photons that travel from r' to r. The number of scattering events per unit energy, per unit time, per unit 
volume at r, along the path 9 is then 
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where bs is the distance between r and r + (Jr, along the path of the photon: 
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Since all angles 9 are equally probable, we compute the angle averaged quantity < S >: 

1 /■'"-,. . ^ S{r',r,E,9) 
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2./n /^T 



<S{r',r,E) >= 
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9max = TT whcn r' < r and 9max = sin^ (^/^') when r' > r. 

In the limit that the wavelength of the incident radiation is much smaller than the Bohr radius, we may assume that 
the radiation sees a gas of electrons and protons, rather than bound atoms. With this approximation, the scattering cross 
section ae-j{Ej) is given by the Klein-Nishina formula (e.g. Peskin & Schroeder 1995) 

a^-^{y)^(TT fa{y), (6) 

where ax is the Thomson cross section for electrons (proton scattering is less important owing to the large mass). f„ is 
given by 
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where y =^ E^ /me- Sabs is the average energy transferred to an electron in a scattering event: 
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Let us make the assumption that the energetic electrons (after scattering) remain in the gas cloud and serve to heat 
up the cloud. As long as low energy photons do not contribute significantly to the heating rate, we may maintain our 
approximation E^ 3> 13.6 eV, while performing the energy integral: 
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H{r) gives the heating rate (energy per unit volume per unit time) at the point r due to primary photons, under the 
approximate assumption that all the primary photon energy is transferred to the gas at the location of the first collision. 
Note, we are not solving the full radiation transfer problem in determining H{r), especially the distributed heating effects 
of secondary pair and photon production, and so our treatment is still an approximate one. 

2.1. The {S{r,r' , E)) function. 
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Fig. 1. — The {S{r, r' , E)) function for m~^ = 100 GeV. S is small for small energies because of the larger scattering cross section. S is also 
small when the point r' is far from r. These effects are more pronounced for small r (high gas density). 



In order to compute the heating rate at a point r using Eq. I10| we need to consider dark matter in a finite volume 
around r. The size of this volume is determined by the function {S{r^ r', E)). When this volume is small compared to the 
size of the halo, one may expect most of the energy released by dark matter annihilation to be absorbed by the gas cloud. 
Figure [T] shows < 5 > for the cases r = 10~^, 10"'^ and 10~^ pc. The three curves in each panel correspond to values of 



r' = 0.999r, r' = 0.99r and r' = 0.9r. As expected, S is largest when r' is close to r. More energetic photons penetrate 
farther through the gas owing to the smaller scattering cross section. 

3. NUMERICAL RESULTS FOR WIMP ANNIHILATION HEATING AND BARYONIC COOLING RATES 

We analyzed the structure of the three highest resolution star- forming minihalos in the numerical simulations presented 
by O'Shea & Norman (2007). We will refer to these simulation runs as A, B, C. These simulations were performed using 
the Enzo code, which is an adaptive mesh refinement (AMR) cosmology code developed by G. Bryan and others (e.g. 
O'Shea et al. 2004). The code couples an N-body particle mesh (PM) solver (Efstathiou et al. 1985; Hockney & Eastwood 
1998) with Eulerian AMR. All simulations were initialized at z = 99 assuming a concordance cosmological model with 
Vim = 0.3, rif, = 0.04, r^cDM = 0.26, r^A = 0.7, h — 0.7, o-s = 0.9 and an Eisenstein & Hu (1999) power spectrum with 
a spectral index n — \. Each of the three star- forming minihalos was formed in a separate simulation with a different 
random seeding of the initial conditions and a box size of 0.3ft.~^ Mpc (comoving). The most massive halo to form in each 
simulation, at z = 15 (typically with a mass ~ lO^M©) was found using a dark matter only computation. The initial 
conditions were then regenerated with both dark matter and baryons such that the Lagrangian volume in which the halo 
formed was resolved at much higher spatial and mass resolution. 

3.1. Dark Matter Density Structure and Heating Rate 

The minimum dark matter particle mass in the simulations was 1.81ft.~^ Mq -^ 2.59 Mq. To avoid effects due to 
the finite size of the dark matter particles, the dark matter density was smoothed on a comoving scale of 0.5 pc, i.e. a 
proper scale of ~ 0.025 pc at z '--^ 20. Note that the gravitational forces due to the baryonic core of the minihalo are not 
smoothed, since there is high spatial resolution in this region. At their final timestep, all simulations had more than 100 
dark matter particles inside a proper scale of 1 pc. 

The dark matter density profile is not resolved in the very central regions of the halo where WIMP annihilation heating 
may be important. Thus the results of these simulations are used as a guide for two extrapolation methods (DM fits #1 
and #2) to estimate the dark matter density on scales of ^ 10"^ — 10^** pc. We fit power laws of the form 

p, = i.(^)"°" (11) 

to the radial dark matter density profiles, evaluated over particular ranges in radius. DM fit #1, shown in Fig. [21 uses only 
data points from the central regions with more than 100 particles per bin (corresponding to scales from about 50 pc down 
to r > 1.5 pc). We find a^ = 1-35, 1.24, 1.65 and A^ = 1674, 1696, 2732 GeV cm'^ for simulations A, B, C, respectively. 
Note that simulation C is has a significantly steeper density profile than A and B. All these density profiles are expected 
to underestimate the dark matter density in the inner regions of the halo, where the baryon density begins to dominate 
over dark matter and adiabatic contraction of the dark matter tends to steepen its profile (Blumenthal et al. 1986). Thus 
we consider DM fit #1 to provide a conservative lower limit to the dark mater density and thus the heating rate. 

Our second extrapolation method, DM fit #2, is based on simulation data from the inner (but not the innermost, 
r < 0.025 pc) regions (see Fig. ^. For simulation A, we have A^ = 2396 GeVcm"'^, a^ = 1.980 when r < 1 pc (and 
A^ = 2409 GeV cm~^, a^ = 1.51 when r > 1 pc). For simulation B, A^ = 1602 GeV cm~^, a^ = 1.946 when r < 1 pc 
(and A^ = 1696 GeVcm"^, a^ = 1.24 when r > 1 pc). For simulation C, we have A^ = 942 GeVcm"^, a-^ = 2.108 
when r < 0.1 pc (and A^ — 2732 GeV cm~'^, a^ = 1.645 when r > 0.1 pc). The density profiles of the inner regions show 
a significant steepening compared to the outer regions. This is likely to be the result of adiabatic contraction, since it 
occurs at a scale '^ 1 pc where the density of baryons begins to dominate. We emphasize again that the simulation points 
on which DM fits #2 are based correspond to bins with only a small (< 100) number of particles and must be treated 
with caution. We do note, however, that on the very innermost scales (~ 0.01 pc) the dark matter density profiles in the 
numerical simulations are even steeper than the analytic DM fits #2. 

If adiabatic contraction of dark matter in the baryon-dominated potential is responsible for sculpting these profiles, 
one expects that the maximum steepness of the dark matter density profile will be equal to that of the baryons, which 
has been shown in simulations to have an approximate power law profile of pgas oc r~^'^ (see i i3.2p . The baryonic density 
profile inevitably flattens in its center, so that the power law profile is only valid inwards to some core radius, Tc. This 
core radius shrinks as collapse proceeds. In situations where we are interested in the global luminosity provided by WIMP 
annihilation in the halo fi J3.3p . we will assume that the dark matter density profile also exhibits this core radius, inside of 
which its density is also constant. 

For the fiducial case, we assume a dark matter particle mass of m^ = 100 GeV, which is typical for a weakly inter- 
acting particle. The photon multiplicity function dN^/dEy for the important annihilation channels is given by the form 
(Bergstrom et al. 1998; Feng et al. 2001) 

dN^ ^ ae^^ ^^2) 

dx x^-^ 

where x — E^/m^ and (a, b) are constants for the particular annihilation channel. We use the values given in Feng, Matchev 
&; Wilczek (2001), and average over the different annihilation channels they considered to obtain (a = 0.9, b = 9.56). We 
set {aav) = 3 X 10~^^ cm'^ s~^ (e.g. Jungman et al. 1996). The heating rate as expressed in equation [TO] was calculated 
for the above density profiles and is shown in Fig. [31 
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Fig. 2. — Dark matter density profiles (spherical averages) for simulations A (top), B (middle), and C (bottom) are shown with the open 
squares. Power law fits to the outer, well-resolved regions (fit #1) are shown with the dashed lines, while fits to the inner (but not innermost 
- see text), Icss-well-resolved, regions (fit #2) are shown with the dotted lines. 



The variation of H{r) on the WIMP mass, m^, is also shown in Fig. [31 We consider cases with m^ = 10 GeV and 1 TeV, 
i.e. factors of 10 below and above our fiducial value. Once the simple m~^ dependence of H{r) (eq. [TU)) is accounted for, 
we see that remaining variations in H{r) are within about a factor of 2. These are due to the m^ dependencies of the 
photon multiplicity function, the energy integral and the {S(r,r' , E)) function. 

It is also informative to compare H{r) to the energy generated by WIMP annihilation per unit volume, per unit time 
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Fig. [3] compares G{r) and H{r) for Simulation A, for DM fit #2. At small distances (r < 100 AU), the two curves are 
very similar, whereas at larger distances, G{r) ^ H{r) since not all the energy generated is absorbed. It is also interesting 
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Fig. 3. — Top panel: WIMP annihilation heating rate, H(r), for simulation A (solid lines), B (long-dashed lines) and C (dot-dashed lines) 
with DM fit #1 (lower, thin lines) and DM fit #2 (upper, thick lines). All these models assume m^ = 100 GeV. The cooling rate, C{r), is 
shown by the dashed line (see i|3.2p . Bottom panel: Variation of H{r) with WIMP mass for simulation A, DM profile #2. WIMP masses 
m^ = 10, 100 and 1000 GeV are shown by the dotted, solid, and dashed lines, respectively. The dot-long-dashed line shows the WIMP 
annihilation energy generation rate, G{r), for m-^ = 100 GeV (it is to be compared with the solid line). At large radii, H(r) becomes much 
smaller than G(r), since not all the energy generated is absorbed. 



to note that H{r) can exceed G{r) at high densities, since unUke G{r), H{r) is not a strictly local function of r. 

3.2. Baryon Density, Temperature, H2 Fraction and Cooling Rate 

In the numerical simulations of O 'Shea & Norman (2007), collapse of baryons was followed with adaptive mesh refinement 
to much higher spatial resolution than the dark matter. Spatial scales down to 115/i~^ — > 164 AU (comoving), i.e. about 
10 AU (proper), were resolved (see Fig. [3]). However, in many cases we consider, knowledge of the gas properties is 
required on even smaller scales. Thus we utilize the results of the even higher resolution smooth particle hydrodynamics 
(SPH) simulations of Yoshida et al. (2006), which were carried out for a single minihalo in which the collapse was followed 
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Fig. 4. — Top panel: radially-averaged H number density, nn, for simulation A (open squares). The solid line shows a power law nn oc r~^'^, 
normalized from the simulations of Yoshida et al. (2006) to have logio(nH/cm~'') = 13.15 at r = 10 AU = 4.88 X 10~^ pc. Middle panel: 
radially-averaged temperature, T, versus radially-averaged H number density, njj for simulation A (open squares). The solid line shows the 
analytic fits estimated from the results of Yoshida et al. (2006). Bottom panel: radially-averaged H2 mass fraction versus radially-averaged 
nil for simulation A (open squares). The solid line shows the analytic fits estimated from the results of Yoshida et al. (2006). 



to higher densities. The structure of baryonic core is controlled mostly by the microphysics of H2 cooling, in which case 
one expects a fairly universal profile of density, temperature and chemical composition. 
For the Hydrogen number density, ripj = 0.76pgas/ 'Tip, we use the analytic fit 

-2.2 



nn = 1-41 X 10' 



UOAuJ 



(14) 



,10AU> 

which is based on conditions at 10 AU in the simulation of Yoshida et al. (2006). Figure [3] shows this fit is also a good 
description of the gas in the simulations of O'Shea & Norman (2007) at r > 200 AU. 

The temperature profile (as a function of density) at densities nn > 10^ cm~^ was fit with an analytic function based 



on the results of Yoshida et al. (2006). For lower densities, our simulations A, B and C are reliable. We use the fitting 
functions (see Fig.g): T ^ 165.5 (nH/cm-3)00754 ^ foj. densities n-a > 10i°cm-3; log(r/i^) = 3-0.025 log^(nH/109cm-3) 
for densities 10^ cm"^ < tih < lO"^" cm"^; T = l78{nii/cm~^)°-°^^ K for densities 10^ cm"^ < n^ < 10^ cm"^; and 
T = 80.5(nH/cni-3)Oi3i5 ^ f^j. densities nn < lO'^ cm-^. 

The molecular hydrogen mass fraction /moi is also fit from the results of the simulations of O'Shea & Norman (2007) 
and Yoshida et al. (2006). We use the fitting functions (see Fig. H]): /moi = 5 x 10~^^(nH/cm~^)°-^^^ for densities 
riH > 10^ cm-3 (Yoshida et al. 2006) and fmo\ = 2.64 x IQ-^inn/cm-^f"^^^ for densities nn < 10^ cm-^ (O'Shea & 
Norman 2007). 

We use the analytic results of HoUenbach & McKee (1979) to implement rotational and vibrational H2 cooling. Collision- 
induced emission (CIE) cooling follows the fit published in Yoshida et al. (2006). We also account for the loss of cooling 
efficiency due to the gas becoming optically thick. The opacity of the gas to cooling radiation was measured from the results 
of Yoshida et al. (2006). The cooling efficiency, 77, follows the fitting function: 77 = 9.4 x 10^(nH/cm~^)~*'-^'*^ for densities 
riH > 10^^ cm~^; rj = 900(nH/cm~^)~°-^^^ for densities 10^^ cm~^ < rin < 10^^ cm~^; 77 = 1.7 x 10^(nH/cm~^)~"-^^^ for 
densities 10^° cm~^ < nn < 10^^ cm"*^; and 77 = 1 for lower densities. 

The cooling rate, C(r), is shown in Fig. [3] for simulation A. The results for simulations B and C are very similar. 

3.3. Comparison of Heating and Cooling Rates and the Equilibrium Structure 

Figure [5] shows the radial profiles of the ratio of the heating rate caused by WIMP annihilation to the gas cooling rate, 
H{r)/C{r), for simulations A, B, C and dark matter density fits #1 and #2, assuming m-^ = 100 GeV. For simulation A, 
the dependence on m^ is also illustrated using m^ = 10 GeV and 1000 GeV. Results are also shown as a function of nn, 
using equation [ni to derive the density at a given radius. 

For the conservative, relatively shallow dark matter density profiles (fit ^1) of simulation A and B and for a fiducial 
WIMP mass of 100 GeV, the heating rate is never larger than ^ 10^^ of the cooling rate for densities up to ?ih — lO^^cm"^. 
In this situation we expect that the star formation process from these gas cores, in particular the protostellar structure 
and accretion rate (e.g. Tan & McKee 2004) would be unaffected by dark matter annihilation heating. The expected mass 
of the stars forming from these halos is then ~ 100 — 2OOM0, set by ionizing feedback, especially disk photoevaporation, 
on the protostellar accretion flow (McKee & Tan 2008). 

The case of DM fit ^1 for simulation C is somewhat steeper (a^ = 1.65) than for simulations A and B (a^ ~ 1.3), and 
this leads to the heating rate becoming comparable to the cooling rate at nn ^ 10^^ cm^'^. This trend is continued for 
the steeper density profiles of DM fit #2, all of which have a central region where WIMP annihilation heating dominates 
over cooling, typically for r < 10^^ pc, i.e. 20 AU, and nn ^ lO""^^ cm~^. One important factor setting this density scale 
is drop in the cooling efficiency of the gas (Yoshida et al. 2006, their Fig. 4) because of the increasing opacity of the gas. 

These results are broadly consistent with the those of Spolyar et al. (2008), showing that the properties of the dark 
matter heating dominated core are relatively insensitive to the details of radiative transport of the WIMP annihilation 
heating that we have included, and that their analytic model for the dark matter density distribution is similar to our 
DM fit #2 to simulated minihalos (see also Freese et al. 2008b). 

The expected consequence of a central core where dark matter heating of the gas dominates its radiative cooling is 
a halt to the collapse. This process would then set the core radius, re, for both the dark matter and baryon density 
distributions. To evaluate this core radius and the total luminosity of the equilibrium structure, we now consider the 
enclosed luminosity profile, L{r), of the halos, defined for WIMP annihilation heating as 

L^(r) = 47r/ dr'r'^ H(r'), (15) 

"'0 

and for baryonic cooling as 

Lgas(r)=47r / dr'r'^Cir'). (16) 

Jo 

The luminosity due to WIMP annihilation will be dominated by the central part of the halo for a^ > 1.5 (eq. [T5)) . which 

is the case for all DM fits #2. We solve for the core radius, r^ for each simulation via the condition L^(r -^ cx)) = 

Lge,sir -^ 00), i.e. 

Airr^ f°° inr^ f°° 

^-H{re) + A-K drr^H{r) = ^C{rc)+in drr^C{r). (17) 

This is an equilibrium condition in the sense that the energy generation rate by dark matter annihilation heating balances 
that radiated away by the baryons. Note that for the profiles analyzed here, both heating and cooling luminosities are 
dominated by the central regions, i.e. the integrals converge to a finite value as r — + 00. Note also that local heating will 
dominate local cooling in an inner region that is somewhat larger than re, and cooling will dominate heating outside of 
this region. 

Figure [H] shows the luminosity profiles for DM fit #2 for simulations A, B and C with m^^ — 100 GeV, and for simulation 
A with TOjj, = 10 GeV and 1 TeV. Also shown is Lgas for the fiducial model. With m^ = 100 GeV, the core radii for 
simulations A, B, C are re = 7.4,2.2,19 AU, respectively, and L^{r — > 00) = 1460, 1540, 781L0, respectively. The dark 
matter density in these cores are p^^c = 1-5 x 10^^,7.5 x 10^^,3.0 x 10^^ GeVcm"'^, respectively. For simulation A 
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Fig. 5. — Ratio of dark matter annihilation heating rate, H{r), to baryonic coohng rate, C{r), for dark matter density fits #1 (top panels) 
and #2 (bottom panels) for simulations A (solid lines), B (long-dashed lines), C (dot-dashed lines) for m~^ = 100 GeV. Results for simulation 
A with m^ = 10 GcV (dotted lines) and m-^ = 1 TcV (dashed lines) arc also shown. 



with m^ = 10 GcV and 1 TeV, we find Tc — 37,0.93 AU, respectively, L^{r -^ cxo) ~ 960, ISIOL©, respectively, and 
p^.c = 6.2 X 10i",9.2 X 10^3 GcVcni^3^ respectively. 

The constancy of these total luminosities is set by the baryonic cooling luminosity, which is dominated by the region 
from ~ 10 AU to ~ 100 AU, and is independent of the much smaller core radius, which adjusts itself so as to provide 
enough WIMP annihilation luminosity to match that radiated away by the baryons. 

We note our solution has made the approximation that the baryonic properties used to estimate Lgas do not include 
the effects of WIMP annihilation heating. This effect would lead to hotter temperatures and more efficient cooling, so 
that the actual equilibrium structure would be somewhat smaller and denser. 

Energy transport from the hotter inner part of the halo to the outer cooler part, either via convection or radiation, 
would be needed to set up a hydrodynamic equilibrium: i.e. a dark matter powered (proto)star. The size of this star 
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Fig. 6. — Enclosed luminosity generated from the heating due to dark matter annihilation, L^ 
(see text). The top panel shows L-^{r) for simulations A (solid line), B (long-dashed line), C (dot-dashed line), assuming m^ 



100 GeV and 



dark matter density fit #2. The same model for simulation A, but with m^ = 10 GeV is shown with the dotted line and with m^ = 1 TeV is 
shown with the dashed line. The bottom panel compares L-^ (r) for simulation A, dark matter density fit #2 and nix ~ ^00 GeV (solid line) 
with the baryonic cooling luminosity, Lgas(»'), (dotted line) for the same halo. 



would be set by the t — 1 surface, beyond which the energy flux from the interior is free to escape. The results of Yoshida 
et al. (2006) imply that the cooling efflciency has dropped to 0.1 by the time nn = 10^^ cm""^, corresponding to about 
30 AU. This is an approximate upper limit to the initial size of the protostar, since the actual gas temperature will be 
somewhat hotter and a larger fraction of the radiated energy would be in the continuum (as opposed to line radiation). 

3.4. Subsequent Evolution of the Dark Matter Powered Protostar 

Since there is a large mass reservoir at large r that is undergoing cooling, baryons will continue to accrete to the central 
protostar, deepening the gravitational potential and thus causing more dark matter to be concentrated here also. For 
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simulation A, DM fit#2 and rri)^ — 100 GeV with Vc = 7.4 AU, the mass inside Tc is 0.17 Mq, and the accretion rate is 
close to O.lAf0yr^^ with infall speeds of about 4 kms^^ (Yoshida et al. 2006), which are mildly supersonic. As the central 
mass grows, the infall is expected to become more supersonic, i.e. closer to the free fall speed of the central protostellar 
mass at its surface. 

The detailed evolution of this dark matter powered protostar is beyond the scope of this paper. Freese et al. (2008c) 
have considered the protostellar structure of such stars, but starting from a more evolved stage when m* = 3 Mq and for 
m* = 2 X 10~^ Mq yr~^. After 4,500 yr of accretion this initial structure has grown in mass to 12 M© and has a radius 
of 2.8 AU and a luminosity of 1.1 x 10^ Lq (ignoring accretion luminosity). The dark matter in the star is assumed to 
increase via adiabatic contraction of the surrounding halo and enough luminosity is provided by dark matter annihilation 
to support stars of at least 1000 Mq, for which L^, ~ 4 x 10® Lq. Being much larger and cooler than protostars on the 
zero age main sequence (ZAMS), these dark matter powered protostars are expected to have relatively weak radiative 
feedback on their accretion envelopes and thus may be able to attain very high masses. 

An alternative possibility is that the baryonic mass of the star grows more rapidly than the enclosed dark matter mass 
and that this more massive star is unable to be supported by WIMP annihilation luminosity. The star contracts and 
reaches central densities and temperatures at which nuclear fusion of H starts. Note that fiducial models of Population 
III.l protostars (e.g. Omukai & Palla 2003; Tan & McKee 2004) do not expect contraction to the main sequence until 
m^ c± 100 Mq, i.e. when the Kelvin-Hclmholz time becomes less than the accretion growth time. The critical feedback 
process, disk photoevaporation, theorized to limit the masses of Pop III.l stars is expected to truncate protostellar growth 
along the ZAMS at about 140 Mq, depending on the accretion rate (McKee & Tan 2008). Thus if WIMP annihilation 
were able to support the protostar only up to ~ 100 Mq, there would be only modest implications for the initial mass 
function of the first stars. 

Understanding the accumulation of dark matter by the protostar is of critical importance for deciding between these 
two scenarios. There are processes and effects that can cause inefficient adiabatic contraction of the dark matter halo in 
response to protostellar accretion. Rapid free-fall collapse of the baryons joining the star-disk system is expected after 
the expansion wave moves out. Standard analytic analyses of adiabatic contraction typically assume very slow changes 
in the baryonic potential (e.g. Blumcnthal et al. 1986). For rapid collapse, one expects a delayed response of the dark 
matter, and during this time jti^dm/'ti* will have decreased. 

Another effect is excitation of the dark matter velocity dispersion by scattering with baryonic density fluctuations. 
These will probably occur most strongly in the accretion disk around the protostar, since such fluctuations are required 
for angular momentum transport, and, if not initially present, will inevitably develop via gravitational instability as the 
disk mass grows (Tan & Blackman 2004). For accretion from a core with a given density and angular momentum, the 
disk size can be described via (Tan & McKee 2004): 

where /kcp is the ratio of the rotational to Keplerian velocities of the gas at the sonic point of the infall and has a fiducial 
value of 0.5 in the simulations of Abel et al. (2002), to*c( is the mass of central protostar and its disk (we typically assume 
a disk mass that is one third of the protostellar mass), and K' is the normalized entropy parameter (related to the density 
of the core at the time of protostar formation: a denser core has a larger value of K'), with fiducial value of 1.0. We see 
that the expected disk sizes when the central mass is 1, 10, 100 Mq are 3.4,66, 1300 AU. Most mass will join the dark 
matter powered protostar via a disk once 771*^ > a few solar masses. 

If the initial protostar fails to accumulate any addition dark matter and retains what it has during Kelvin-Helmholz 
contraction (which is expected since the contraction time is much longer than the dynamical time), then, ignoring deple- 
tion (discussed below, but which occurs on timescales much longer than the initial accretion growth time), the WIMP 
annihilation luminosity of the star is approximately described by 

Lx^L^.o(^) , (19) 

where L^^o — 10'^ Lq and r*.o — 10 AU are fiducial values for the initial luminosity and radius of the star, respectively, 
and assuming homologous evolution of the density profile. The protostellar size in models with X^ = is expected to be 
~ 1 AU for m* < 10 Mq (Tan & McKee 2004). Compressing a protostar to this size would lead to L^ 2± 10® Lq for the 
above starting conditions, which is enough luminosity to support a ~ 100 Mq star. So in this case the protostar would 
be able to be supported by WIMP annihilation giving it > AU sizes until about 100 Mq . The ZAMS radius of a 100 Mq 
star is about 4 Rq = 0.02 AU, so it is clear that contraction to this scale would not be possible until very large stellar 
masses, unless the WIMP content of the protostar became depleted. 

These conclusions are of course sensitive to the initial condition. We expect L^^ ~ 10"^ L© to be relatively robust as it 
is set by the baryonic cooling properties of the halo, but the initial size of the dark matter core that provides this WIMP 
annihilation luminosity depends quite sensitively on the uncertain WIMP mass and the spatial distribution of the dark 
matter. For example, if the dark matter core has an initial size of about 1 AU (e.g. similar to simulation B, DM fit#2, 
m-^ = 100 GeV and models with m^ = 1 TcV) and assuming the initial protostellar size scales proportionately and is 
also about equal to this size, then the potential luminosity available from contraction of the protostar to a given radius 
would be a thousand times less than the case considered above. Nevertheless, a contraction to 10 Rq from 1 AU, would 
still lead to a WIMP annihilation luminosity of 8 x 10® Lq. 
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We now compare the growth time of the protostar (i.e. the time since its formation, its age), 
U = 2.92 X 104ii:'-i5/7(m,/100Mo)io/7 yr (Tan & McKee 2004) to the WIMP depletion time: 

-1 



^dcp . 



Px 
Px 



105 



Px 



100 GeV VlQi^GeVcm- 



Myr. 



(20) 



Px < '^aV > 

If the protostar contracts from an initial radius of 10 AU (where we find p^ ~ 10^^ GeV cm^'^, see above) to 1 AU then 
idep — 10^ yr. We see that, if replenishment of WIMPs in the protostar is negligible, then depletion can become important 
for AU scale protostars of ~ 100 Mq. These timescales are sensitive to m^ and the dark matter density in the initial 
equilibrium core, which shows significant variation between simulation A, B, and C. 



4. CONCLUSIONS 

We have investigated the effects of WIMP dark matter annihilation on the formation of Population III.l stars by 
analyzing the results of cosmological simulations that follow the gravitational collapse of baryons and dark matter. While 
these simulations (O'Shea & Norman 2007; Yoshida et al. 2006) have followed the baryons to very high densities at 
scales < 1 AU, the dark matter is only well-resolved down to scales ~ 1 pc. Thus we have considered various power law 
extrapolations of the dark matter density profile towards the center. 

If one assumes the dark matter profile is self-similar, extending inwards from the dark matter dominated regime with 
p^ ex r~°'^ and a^ ~ 1.5, then, for a fiducial WIMP mass of 100 GeV, the dark matter annihilation heating is typically 
negligible {^ 10~ of the cooling rate). This conclusion would be unchanged for a reduction in the WIMP mass by a 
factor of 10 or more. One of the simulated minihalos (G) exhibits a slightly steeper density profile in the well-resolved 
region (a^ = 1.65, and in this case WIMP annihilation heating does become important at nn ^ 10^^ cm^^. 

However, there are theoretical reasons to expect that the dark matter density profile will steepen because of adiabatic 
contraction in the baryon-dominated core. Indeed, this process appears to be occurring in the simulations of O'Shea & 
Norman (2007), although it is not well-resolved. A value of a^ ~ 2.0 appears to be a better description of the dark matter 
density profile in this region. For such a profile, and again for a 100 GeV particle, the dark matter annihilation heating 
now exceeds baryonic cooling for densities nn > lO-'^'^cm"'^, in agreement with the previous study of Spolyar et al. (2008). 

We considered the properties of equilibrium halos in which the density distributions of the baryons and dark matter 
exhibit a constant density central core. The luminosity that is generated is ~ 10'^ i© and is relatively invariant, being set 
by baryonic cooling processes. The sizes of the central cores range from ^^ 1 to 40 AU. 

The detailed effects of this extra heating on the protostellar structure remain to be determined. We expect that 
subsequent baryonic growth of the protostar will occur more rapidly than its accumulation of dark matter, because the 
baryons are undergoing rapid free fall collapse followed by disk accretion. However, even if the initial protostar does not 
gain any additional dark matter, its initial dark matter content in the subsolar mass core could be sufficient to prevent 
contraction to the zero age main sequence for masses of 100 M© or greater. Such circumstances could have dramatic 
implications for the masses of Pop III.l stars, conceivably raising the mass scale to a regime important for the formation 
of supermassive black holes. These conclusions depend sensitively on the initial protostellar core, which now needs to 
be studied with self-consistent cosmological simulations that include the influence of WIMP annihilation heating on the 
baryons. 

We thank T. Abel, K. Freese, G. McKee, D. Spolyar, and N. Yoshida for discussions. The research of JGT is sup- 
ported by NSF CAREER grant AST-0645412. A.N. acknowledges financial support from the Deutsche Forschungsge- 
meinschaft(DFG) International Research Training Group GRK 881. 

REFERENCES 



Abel, T., Bryan, G. L, & Norman, M. L. 2002, Science, 295, 93 

Ascasibar, Y., 2007, A&A, 462, L65 

Bergstrom, L., Ullio, P. & Buckley, J. H., 1998, Astroparticle Physics, 

9, 137 
Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, 

ApJ, 301, 27 
Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23 
Efstathiou, G., Davis, M., White, S. D. M., & Frenk, C. S., 1985, 

ApJS, 57, 241 
Eisenstein, D. J. & Hu, W., 1999, ApJ, 511, 5 
Feng, J. L., Matchev, K. T., & Wilczek, F., 2001, Phys. Rev. D63, 

045024 
Freese, K., Bodelheimer, P., Spolyar, D. & Gondolo, P. 2008c, arXiv: 

0806.0617 
Freese, K., Gondolo, P., Sellwood, J. A. & Spolyar, D. 2008b, arXiv: 

0805.3540 
Freese, K., Spolyar, D. & Aguirre, A. 2008a, arXiv: 0802.1724 
Hall, J. & Gondolo, P., 2006, Phys. Rev. D74, 063511 
Hockney, R. W. & Eastwood, J. W., 1998, Computer Simulation 

using Particles (Bristol: lOP) 
HoUenbach, D. & McKee, C. F. 1979, ApJS, 41, 555 
locco, F., 2008, ApJ, 677, LI 
Jungman, G., Kamionkowski, M., & Griest, K., 1996, Physics 

Reports, 267, 195 



Mckee, C. F., & Tan, J. C, 2008, ApJ, 681, 771 
Navarro, J. F., Frenk, C. S. & White, S. 1996, ApJ, 462, 563 
Omukai, K., & Palla, F. 2003, ApJ, 589, 677 

O'Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., 
Harkness, R., & Kritsuk, A. 2004, in "Adaptive Mesh Refinement - 
Theory and Applications", Eds. T. Plewa, T. Linde & V. G. Weirs, 
Springer Lecture Notes in Computational Science and Engineering, 
2004 ( astro-ph/0403044 ) 
O'Shea, a. W., McKee, C. F., Heger, A., & Abel, T. 2008, in First 
Stars III, eds. B. W. O'Shea, A. Heger & T. Abel (ar Xiv:0801.2124ll 
O'Shea, B. W. & Norman, M. L., 2007, ApJ, 654, 66 
Peskin, M. & Schroeder, D. V., 1995, Westview Press 
Ripamonti, E., Mapelli, M., & Ferrara, A., 2007, MNRAS, 375, 1399 
Spolyar, D., Freese, K., & Gondolo, P., 2008, Physical Review 

Letters, 100, 051101 
Tan, J. C, & Blackman, E. G. 2004, ApJ, 603, 401 
Tan, J. C, & Mckee, C. F. 2004, ApJ, 603, 383 



Taoso, M., Bertone, G., Meynet, G. 

arXiv:0806.2681 
Yoshida, N., Omukai, K., Hernquist, L., & 

6 
Yoon , S.-C, lo cco, F., Akiyama, S. 

iarXiv:0806.2662i 



Ekstrom, S. 2008, 
Abel, T., 2006, ApJ, 652, 
2008, ApJ, submitted. 



